####Model results####
##a) alternative aggregation that mirrors RAI self-rule aggregation (table X24, models A80-A83)
r3_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_rai","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r3_m1g_t1_sc_cse <- data.frame(cluster.se(r3_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r3_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_rai","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r3_m1g_t1_sf_cse <- data.frame(cluster.se(r3_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r3_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_rai","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3_m1d_t1_sb_cse <- data.frame(cluster.se(r3_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r3_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_rai*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3_m2d_t1_sb_cse <- data.frame(cluster.se(r3_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r3_m1g_t1_sc, r3_m1g_t1_sf, r3_m1d_t1_sb, r3_m2d_t1_sb, se=c(r3_m1g_t1_sc_cse, r3_m1g_t1_sf_cse, r3_m1d_t1_sb_cse, r3_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A80)", "Min. (model A81","Maj./min. dyad (model A82)","Maj./min. dyad (model A83)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X24. Additional results: alternative aggregation of territorial autonomy index (RAI-style additive index).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex24.txt")
stargazer(r3_m1g_t1_sc, r3_m1g_t1_sf, r3_m1d_t1_sb, r3_m2d_t1_sb, se=c(r3_m1g_t1_sc_cse, r3_m1g_t1_sf_cse, r3_m1d_t1_sb_cse, r3_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A80)", "Min. (model A81","Maj./min. dyad (model A82)","Maj./min. dyad (model A83)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X24. Additional results: alternative aggregation of territorial autonomy index (RAI-style additive index).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex24.html")
##b) every component required (table X25, models A84-A87)
r2b_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_df2","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r2b_m1g_t1_sc_cse <- data.frame(cluster.se(r2b_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r2b_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_df2","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r2b_m1g_t1_sf_cse <- data.frame(cluster.se(r2b_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r2b_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_df2","included_excluded", "excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r2b_m1d_t1_sb_cse <- data.frame(cluster.se(r2b_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r2b_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_df2*included_excluded", "excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r2b_m2d_t1_sb_cse <- data.frame(cluster.se(r2b_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r2b_m1g_t1_sc, r2b_m1g_t1_sf, r2b_m1d_t1_sb, r2b_m2d_t1_sb, se=c(r2b_m1g_t1_sc_cse, r2b_m1g_t1_sf_cse, r2b_m1d_t1_sb_cse, r2b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A84)", "Min. (model A85)","Maj./min. dyad (model A86)","Maj./min. dyad (model A87)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X25. Additional results: alternative aggregation of territorial autonomy index (all components required).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex25.txt")
stargazer(r2b_m1g_t1_sc, r2b_m1g_t1_sf, r2b_m1d_t1_sb, r2b_m2d_t1_sb, se=c(r2b_m1g_t1_sc_cse, r2b_m1g_t1_sf_cse, r2b_m1d_t1_sb_cse, r2b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A84)", "Min. (model A85)","Maj./min. dyad (model A86)","Maj./min. dyad (model A87)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X25. Additional results: alternative aggregation of territorial autonomy index (all components required).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex25.html")
##c) fiscal autonomy required (table X26, models A88-A91)
r2a_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_df","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r2a_m1g_t1_sc_cse <- data.frame(cluster.se(r2a_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r2a_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_df","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r2a_m1g_t1_sf_cse <- data.frame(cluster.se(r2a_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r2a_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_df","included_excluded", "excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r2a_m1d_t1_sb_cse <- data.frame(cluster.se(r2a_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r2a_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_df*included_excluded", "excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r2a_m2d_t1_sb_cse <- data.frame(cluster.se(r2a_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r2a_m1g_t1_sc, r2a_m1g_t1_sf, r2a_m1d_t1_sb, r2a_m2d_t1_sb, se=c(r2a_m1g_t1_sc_cse, r2a_m1g_t1_sf_cse, r2a_m1d_t1_sb_cse, r2a_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A88)", "Min. (model A89)","Maj./min. dyad (model A90)","Maj./min. dyad (model A91)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X26. Additional results: alternative aggregation of territorial autonomy index (fiscal autonomy required).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex26.txt")
stargazer(r2a_m1g_t1_sc, r2a_m1g_t1_sf, r2a_m1d_t1_sb, r2a_m2d_t1_sb, se=c(r2a_m1g_t1_sc_cse, r2a_m1g_t1_sf_cse, r2a_m1d_t1_sb_cse, r2a_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A88)", "Min. (model A89)","Maj./min. dyad (model A90)","Maj./min. dyad (model A91)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X26. Additional results: alternative aggregation of territorial autonomy index (fiscal autonomy required).", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex26.html")
##d) special vs statewide autonomy (table X27, models A92-A95)
r3g_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_special","sa_territory_t_swide","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r3g_m1g_t1_sc_cse <- data.frame(cluster.se(r3g_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r3g_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_territory_t_special","sa_territory_t_swide","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r3g_m1g_t1_sf_cse <- data.frame(cluster.se(r3g_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r3g_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_special","sa_territory_t_swide","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3g_m1d_t1_sb_cse <- data.frame(cluster.se(r3g_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r3g_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_territory_t_special*included_excluded","sa_territory_t_swide*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3g_m2d_t1_sb_cse <- data.frame(cluster.se(r3g_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r3g_m1g_t1_sc, r3g_m1g_t1_sf, r3g_m1d_t1_sb, r3g_m2d_t1_sb, se=c(r3g_m1g_t1_sc_cse, r3g_m1g_t1_sf_cse, r3g_m1d_t1_sb_cse, r3g_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A92)", "Min. (model A93)","Maj./min. dyad (model A94)","Maj./min. dyad (model A95)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X27. Additional results: special autonomy arrangements vs. state-wide arrangements.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (special)","Territorial autonomy (state-wide)","Territorial autonomy (special) x included/excluded","Territorial autonomy (state-wide) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex27.txt")
stargazer(r3g_m1g_t1_sc, r3g_m1g_t1_sf, r3g_m1d_t1_sb, r3g_m2d_t1_sb, se=c(r3g_m1g_t1_sc_cse, r3g_m1g_t1_sf_cse, r3g_m1d_t1_sb_cse, r3g_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A92)", "Min. (model A93)","Maj./min. dyad (model A94)","Maj./min. dyad (model A95)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X27. Additional results: special autonomy arrangements vs. state-wide arrangements.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Territorial autonomy (special)","Territorial autonomy (state-wide)","Territorial autonomy (special) x included/excluded","Territorial autonomy (state-wide) x included/excluded","Included","Included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex27.html")

####Figure A15####
##a) alternative aggregation that mirrors RAI self-rule aggregation (civil)
me_r3_m1g_t1_sc <- margins_summary(r3_m1g_t1_sc, variables = "sa_territory_t_rai", vcov = cluster.vcov(r3_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r3_m1g_t1_sc$term <- "second-order maj."
me_r3_m1g_t1_sf <- margins_summary(r3_m1g_t1_sf, variables = "sa_territory_t_rai", vcov = cluster.vcov(r3_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r3_m1g_t1_sf$term <- "second-order min."
me_r3_m1g <- rbind.fill(me_r3_m1g_t1_sc, me_r3_m1g_t1_sf)
me_r3_m1g$term <- as.factor(me_r3_m1g$term)
me_r3_m1g$model <- "alternative aggregation"
me_r3_m1g$estimate <- me_r3_m1g$AME
me_r3_m1g$conf.low <- me_r3_m1g$lower
me_r3_m1g$conf.high <- me_r3_m1g$upper
##b) alternative aggregation that mirrors RAI self-rule aggregation (communal)
me_r3_m1d_t1_sb <- margins_summary(r3_m1d_t1_sb, variables = "sa_territory_t_rai", vcov = cluster.vcov(r3_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3_m1d_t1_sb$term <- "all"
me_r3_m1d_t1_sb$model <- "alternative aggregation"
me_r3_m2d_t1_sb <- margins_summary(r3_m2d_t1_sb, variables = "sa_territory_t_rai", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r3_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3_m2d_t1_sb$term <- ifelse(me_r3_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r3_m2d_t1_sb$model <- "alternative aggregation"
me_r3_m13d <- rbind.fill(me_r3_m1d_t1_sb, me_r3_m2d_t1_sb)
me_r3_m13d <- subset(me_r3_m13d, term != "other")
me_r3_m13d$term <- as.factor(me_r3_m13d$term)
me_r3_m13d$estimate <- me_r3_m13d$AME
me_r3_m13d$conf.low <- me_r3_m13d$lower
me_r3_m13d$conf.high <- me_r3_m13d$upper
##c) every component required (civil)
me_r2b_m1g_t1_sc <- margins_summary(r2b_m1g_t1_sc, variables = "sa_territory_t_df2", vcov = cluster.vcov(r2b_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r2b_m1g_t1_sc$term <- "second-order maj."
me_r2b_m1g_t1_sf <- margins_summary(r2b_m1g_t1_sf, variables = "sa_territory_t_df2", vcov = cluster.vcov(r2b_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r2b_m1g_t1_sf$term <- "second-order min."
me_r2b_m1g <- rbind.fill(me_r2b_m1g_t1_sc, me_r2b_m1g_t1_sf)
me_r2b_m1g$term <- as.factor(me_r2b_m1g$term)
me_r2b_m1g$model <- "all autonomy components required"
me_r2b_m1g$estimate <- me_r2b_m1g$AME
me_r2b_m1g$conf.low <- me_r2b_m1g$lower
me_r2b_m1g$conf.high <- me_r2b_m1g$upper
##d) every component required (communal)
me_r2b_m1d_t1_sb <- margins_summary(r2b_m1d_t1_sb, variables = "sa_territory_t_df2", vcov = cluster.vcov(r2b_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r2b_m1d_t1_sb$term <- "all"
me_r2b_m1d_t1_sb$model <- "all autonomy components required"
me_r2b_m2d_t1_sb <- margins_summary(r2b_m2d_t1_sb, variables = "sa_territory_t_df2", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r2b_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r2b_m2d_t1_sb$term <- ifelse(me_r2b_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r2b_m2d_t1_sb$model <- "all autonomy components required"
me_r2b_m13d <- rbind.fill(me_r2b_m1d_t1_sb, me_r2b_m2d_t1_sb)
me_r2b_m13d <- subset(me_r2b_m13d, term != "other")
me_r2b_m13d$term <- as.factor(me_r2b_m13d$term)
me_r2b_m13d$estimate <- me_r2b_m13d$AME
me_r2b_m13d$conf.low <- me_r2b_m13d$lower
me_r2b_m13d$conf.high <- me_r2b_m13d$upper
##e) fiscal autonomy required (civil)
me_r2a_m1g_t1_sc <- margins_summary(r2a_m1g_t1_sc, variables = "sa_territory_t_df", vcov = cluster.vcov(r2a_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r2a_m1g_t1_sc$term <- "second-order maj."
me_r2a_m1g_t1_sf <- margins_summary(r2a_m1g_t1_sf, variables = "sa_territory_t_df", vcov = cluster.vcov(r2a_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r2a_m1g_t1_sf$term <- "second-order min."
me_r2a_m1g <- rbind.fill(me_r2a_m1g_t1_sc, me_r2a_m1g_t1_sf)
me_r2a_m1g$term <- as.factor(me_r2a_m1g$term)
me_r2a_m1g$model <- "fiscal autonomy required"
me_r2a_m1g$estimate <- me_r2a_m1g$AME
me_r2a_m1g$conf.low <- me_r2a_m1g$lower
me_r2a_m1g$conf.high <- me_r2a_m1g$upper
##f) fiscal autonomy required (communal)
me_r2a_m1d_t1_sb <- margins_summary(r2a_m1d_t1_sb, variables = "sa_territory_t_df", vcov = cluster.vcov(r2a_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r2a_m1d_t1_sb$term <- "all"
me_r2a_m1d_t1_sb$model <- "fiscal autonomy required"
me_r2a_m2d_t1_sb <- margins_summary(r2a_m2d_t1_sb, variables = "sa_territory_t_df", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r2a_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r2a_m2d_t1_sb$term <- ifelse(me_r2a_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r2a_m2d_t1_sb$model <- "fiscal autonomy required"
me_r2a_m13d <- rbind.fill(me_r2a_m1d_t1_sb, me_r2a_m2d_t1_sb)
me_r2a_m13d <- subset(me_r2a_m13d, term != "other")
me_r2a_m13d$term <- as.factor(me_r2a_m13d$term)
me_r2a_m13d$estimate <- me_r2a_m13d$AME
me_r2a_m13d$conf.low <- me_r2a_m13d$lower
me_r2a_m13d$conf.high <- me_r2a_m13d$upper
##g) special vs statewide autonomy (civil)
me_r3g_m1g_t1_sc <- margins_summary(r3g_m1g_t1_sc, variables = c("sa_territory_t_special","sa_territory_t_swide"), vcov = cluster.vcov(r3g_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r3g_m1g_t1_sc$term <- "second-order maj."
me_r3g_m1g_t1_sf <- margins_summary(r3g_m1g_t1_sf, variables = c("sa_territory_t_special","sa_territory_t_swide"), vcov = cluster.vcov(r3g_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r3g_m1g_t1_sf$term <- "second-order min."
me_r3g_m1g <- rbind.fill(me_r3g_m1g_t1_sc, me_r3g_m1g_t1_sf)
me_r3g_m1g$term <- as.factor(me_r3g_m1g$term)
me_r3g_m1g$model <- ifelse(me_r3g_m1g$factor == "sa_territory_t_special", "special autonomy", "statewide arrangement")
me_r3g_m1g$estimate <- me_r3g_m1g$AME
me_r3g_m1g$conf.low <- me_r3g_m1g$lower
me_r3g_m1g$conf.high <- me_r3g_m1g$upper
##h) special vs statewide autonomy (communal)
me_r3g_m1d_t1_sb <- margins_summary(r3g_m1d_t1_sb, variables = c("sa_territory_t_special","sa_territory_t_swide"), vcov = cluster.vcov(r3g_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3g_m1d_t1_sb$term <- "all"
me_r3g_m1d_t1_sb$model <- ifelse(me_r3g_m1d_t1_sb$factor == "sa_territory_t_special", "special autonomy", "statewide arrangement")
me_r3g_m2d_t1_sb <- margins_summary(r3g_m2d_t1_sb, variables = c("sa_territory_t_special","sa_territory_t_swide"), at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r3g_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3g_m2d_t1_sb$term <- ifelse(me_r3g_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r3g_m2d_t1_sb$model <- ifelse(me_r3g_m2d_t1_sb$factor == "sa_territory_t_special", "special autonomy", "statewide arrangement")
me_r3g_m13d <- rbind.fill(me_r3g_m1d_t1_sb, me_r3g_m2d_t1_sb)
me_r3g_m13d <- subset(me_r3g_m13d, term != "other")
me_r3g_m13d$term <- as.factor(me_r3g_m13d$term)
me_r3g_m13d$estimate <- me_r3g_m13d$AME
me_r3g_m13d$conf.low <- me_r3g_m13d$lower
me_r3g_m13d$conf.high <- me_r3g_m13d$upper
##i) combined (civil)
figurea15_plot_data_g <- rbind.fill(me_r3_m1g,me_r2a_m1g, me_r2b_m1g,me_r3g_m1g)
figurea15_plot_data_g$model <- as.factor(figurea15_plot_data_g$model)
figurea15_plot_data_g$model = factor(figurea15_plot_data_g$model,levels(figurea15_plot_data_g$model)[c(2,1,3,5,4)])
figurea15_plot_data_g <- figurea15_plot_data_g[order(figurea15_plot_data_g$model),]
figurea15_plot_g <- dwplot(figurea15_plot_data_g,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("a) civil violence") +
  scale_color_manual(name = "Coefficient for:", values = c("#7fc97f","#7fc97f","#beaed4","#fdc086","#386cb0"))+ scale_linetype_manual(name = "Coefficient for:", values = c(5,4,3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(21,16,17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of civil violence") + ylab("group type")+ #+
  guides(shape = guide_legend(nrow=3,"model",reverse=TRUE), colour = guide_legend(nrow=3,"model",reverse=TRUE),linetype = guide_legend(nrow=3,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format())
##i) combined (communal)
figurea15_plot_data_d <- rbind.fill(me_r3g_m13d,me_r2a_m13d,me_r2b_m13d,me_r3_m13d)
figurea15_plot_data_d$model <- as.factor(figurea15_plot_data_d$model)
figurea15_plot_data_d$model = factor(figurea15_plot_data_d$model,levels(figurea15_plot_data_d$model)[c(2,1,3,5,4)])
figurea15_plot_data_d <- figurea15_plot_data_d[order(figurea15_plot_data_d$model),]
figurea15_plot_d <- dwplot(figurea15_plot_data_d,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("b) communal violence maj./min. dyad") +
  scale_color_manual(name = "Coefficient for:", values = c("#7fc97f","#7fc97f","#beaed4","#fdc086","#386cb0"))+ scale_linetype_manual(name = "Coefficient for:", values = c(5,4,3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(21,16,17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type")+ #+
  guides(shape = guide_legend(nrow=3,"model",reverse=TRUE), colour = guide_legend(nrow=3,"model",reverse=TRUE),linetype = guide_legend(nrow=3,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format(), breaks = c(0, 0.0025, 0.005))
##k) combined (civil)
figurea15 <- grid_arrange_shared_legend(figurea15_plot_g, figurea15_plot_d, ncol=2, nrow=1)
ggsave(figurea15, file='../figures/figurea15.pdf', width = 14, height = 6.5, units="cm",dpi=1000)